%%
 
figure(1)
close
figure(1)
hold on

plot_OPAMS_RMSD='off';

%Ternaries2plot = [4 3 1; 2 3 1; 4 2 3];
Ternaries2plot = [2 3 1;  4 3 1; 4 2 3];
%Ternaries2plot = [4 3 1];
%Ternaries2plot = [4 2 3];
%Ternaries2plot = [2 3 1];
Verticies = {'Qtz' 'Plag' 'Oliv' 'Cpx'};
NamesVerticies = {'QTZ' 'PLAG' 'OLIVINE' 'CPX'};


t1 = Ternary2XY([.3 .7 0 ]);
axishere = [...
    t1(1) 1 t1(2) 0.866
    0 1 0 0.866
    0 1 0 0.866];



axishere = [...
    0 1 0 0.866
    0 1 0 0.866
    0 1 0 0.866];



i=1;
for numTernaryPlots = 1:size(Ternaries2plot,1)
    IndiciesHere = Ternaries2plot(numTernaryPlots,:);
    NameHere = [add2NAME,'_',NamesVerticies{IndiciesHere(1)},'-',NamesVerticies{IndiciesHere(2)},'-',NamesVerticies{IndiciesHere(3)}];
    

    subaxis(1,3,numTernaryPlots,'Margin',.02,'Spacing',.05)
    hold on
    
    
    snippet =  ['''ko'',''linewidth'',1,''MarkerSize'',1,''MarkerFaceColor'',grey7'',''Color'',grey7'''];
    GenericTernaryPlotter5(Galeetal2013_Glasses_Tormey, IndiciesHere, Verticies,snippet,[],[],'no');
    
    
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
        % plots nearest MSP of erupted basalts

    %      for ij =highlight_ReverseFC; %1:size(allMSPs_erupted,3)
    %         GenericTernaryPlotter4(allMSPs_erupted(:,IndiciesHere(1)+2,ij),allMSPs_erupted(:,IndiciesHere(2)+2,ij),allMSPs_erupted(:,IndiciesHere(3)+2,ij), Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
    %             '''.'',''linewidth'',1,''Color'',azure',[],[],'no')
    %         labels2 = cellstr(num2str([(round(allMSPs_erupted(:,1,ij),1))],2));
    % %         if ij ==1
    % %         GenericTernaryPlotter4(allMSPs_erupted(:,IndiciesHere(1)+2,ij),allMSPs_erupted(:,IndiciesHere(2)+2,ij),allMSPs_erupted(:,IndiciesHere(3)+2,ij), Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
    % %             '''.'',''linewidth'',1,''Color'',azure',[],[],'no',labels2(allMSPs_transitionIndicies(ij,:)),allMSPs_transitionIndicies(ij,:),'''FontWeight'',''bold'',''FontSize'',15,''FontAngle'',''italic''')
    % %         end
    %      end
    

    
    % PLOTs other saturation boundaries:

    %         GenericTernaryPlotter4(OPAM(highlight_ReverseFC,IndiciesHere(1)),OPAM(highlight_ReverseFC,IndiciesHere(2)),OPAM(highlight_ReverseFC,IndiciesHere(3)), Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
    %             '''d'',''linewidth'',1,''Color'',black,''MarkerFaceColor'',black,''MarkerSize'',eval(runningSizesBad{1}),''HandleVisibility'',''on''',[],[],'no')
    %
    %
    %         GenericTernaryPlotter4(OPAMout(:,IndiciesHere(1)),OPAMout(:,IndiciesHere(2)),OPAMout(:,IndiciesHere(3)), Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
    %             '''.'',''linewidth'',1,''Color'',black,''MarkerFaceColor'',black,''MarkerSize'',eval(runningSizesBad{1}),''HandleVisibility'',''on''',[],[],'no')
    %
    %
    %         GenericTernaryPlotter4(OPALMout(:,IndiciesHere(1)),OPALMout(:,IndiciesHere(2)),OPALMout(:,IndiciesHere(3)), Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
    %             '''.'',''linewidth'',1,''Color'',green,''MarkerFaceColor'',black,''MarkerSize'',eval(runningSizesBad{1}),''HandleVisibility'',''on''',[],[],'no')
    %
    
    
    %             for ij =highlight_ReverseFC
    %         GenericTernaryPlotter4(OPALM(:,IndiciesHere(1)+2,ij),OPALM(:,IndiciesHere(2)+2,ij),OPALM(:,IndiciesHere(3)+2,ij), Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
    %             '''h'',''markersize'',15,''Color'',red,''HandleVisibility'',''off''',[],[],'no')
    %             end
    %
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   if strcmp(plot_OPAMS_RMSD ,'on') == 1
    % plots MSPs
            GenericTernaryPlotter4(squeeze(allMSPs(:,IndiciesHere(1)+2,highlight_ReverseFC)),...
           squeeze(allMSPs(:,IndiciesHere(2)+2,highlight_ReverseFC)),squeeze(allMSPs(:,IndiciesHere(3)+2,highlight_ReverseFC)),...
           Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
            '''.'',''linewidth'',1,''Color'',black,''HandleVisibility'',''off''',[],[],'no')
    % The line error between corrected and true primary melts
    zzz=1;
    for k = highlight_ReverseFC
        %DataAndModel(1,:)= runningTormey(subdata(k),:);
        DataAndModel(1,:) = dataX(k,:);
        DataAndModel(2,:) = modelX(k,:);
        labels = {};
        %snippet = [ ''':'',''linewidth'',3,''MarkerSize'',8,''Color'',',ColorsHere_all{subdata(k)},''',''HandleVisibility'',''off'''];
        if zzz==1
            snippet =  ['''c-'',''linewidth'',1,''MarkerSize'',1,''HandleVisibility'',''on'''];
        else
            snippet =  ['''c-'',''linewidth'',1,''MarkerSize'',1,''HandleVisibility'',''off'''];
        end
        
        GenericTernaryPlotter4(DataAndModel(:,IndiciesHere(1)),DataAndModel(:,IndiciesHere(2)),DataAndModel(:,IndiciesHere(3)), Verticies{IndiciesHere(1)}, Verticies{IndiciesHere(2)},Verticies{IndiciesHere(3)},...
            snippet ,[],[],'no',labels,1,'''FontWeight'',''bold'',''FontSize'',11,''HorizontalAlignment'',''center'',''VerticalAlignment'',''bottom''')
        %  '''o-k'',''linewidth'',1,''MarkerSize'',8',[],[],'no',labels{k},1,'''FontWeight'',''bold'',''FontSize'',11,''HorizontalAlignment'',''center'',''VerticalAlignment'',''bottom''')
        zzz=zzz+1;
    end
    aaas=1;
    for zz = UniqueLabels_Rows'
        i =find(ismember(CombinedLabel, CombinedLabel(zz)));         
        
        p=ismember(i,highlight_ReverseFC);
        i(p==0)=[];

        % The closest true primary melts
        snippet =  ['''',runningMarkers{zz},''',''linewidth'',1.5,''MarkerSize'',8,''MarkerFaceColor'',''cyan'',''MarkerEdgeColor'',', runningColors_FILL{llkm},',''HandleVisibility'',''off'''];
        GenericTernaryPlotter5(modelX(i,:), IndiciesHere, Verticies,snippet,[],[],'no');
     
        aaas=aaas+1;
    end
   end
 %     plot(0,0,'o','Color','none')
%     plot(0,0,'o','Color','none')
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
        
    
    



    % FC PATH Have to plot each of these as its own line
    a=1;
    for llkm=highlight_ReverseFC
        if a==1
        snippet =  ['''r-'',''linewidth'',3,''MarkerSize'',1,''Color'',',runningColors_FILL{llkm}];
        else
        snippet =  ['''r-'',''linewidth'',3,''MarkerSize'',1,''Color'',',runningColors_FILL{llkm},''',''HandleVisibility'',''off'''];
        end
        GenericTernaryPlotter5(BestFit_FCPath_Tormey_n{llkm}, IndiciesHere, Verticies,snippet,[],[],'no');
        a=a+1;
    end
    
 
    % Original Data that RevPet does not find a Valid Solution
    aaas=1;
    for zztemp = 1:size(Unique_BAD_labels_Rows,1)    
        Unique_BAD_label_here = Unique_BAD_labels{zztemp};
        i =find(ismember(CombinedLabel, Unique_BAD_label_here));
        p=ismember(i,BAD_ReverseFC);
        i(p==0)=[];
        
        zz = i(1);

      %  snippet=['''',runningMarkers{zz},''',''linewidth'',1,''MarkerSize'',',runningSizesBad{1},',''MarkerEdgeColor'',',runningColors_LINE{zz},''',''MarkerFaceColor'',',runningColors_FILL{zz},''',''HandleVisibility'',''on'''];
        snippet =  ['''',runningMarkers{zz},''',''linewidth'',1,''MarkerSize'',',runningSizesBad{1},',''MarkerEdgeColor'',', runningColors_LINE{zz},',''MarkerFaceColor'',', runningColors_FILL{zz},''',''HandleVisibility'',''on'''];
        
        GenericTernaryPlotter5(runningTormey((i),:), IndiciesHere, Verticies,snippet,[],[],'no');
        aaas=aaas+1;
    end
    
 
   % Original Data with Valid Solutions
    aaas=1;
    for zz = UniqueLabels_Rows'
        i =find(ismember(CombinedLabel, CombinedLabel(zz)));
        p=ismember(i,highlight_ReverseFC);
        i(p==0)=[];
        
        snippet =  ['''',runningMarkers{zz},''',''linewidth'',1,''MarkerSize'',',runningSizes{zz},',''MarkerEdgeColor'',', runningColors_LINE{zz},',''MarkerFaceColor'',', runningColors_FILL{zz},''',''HandleVisibility'',''on'''];
        GenericTernaryPlotter5(runningTormey((i),:), IndiciesHere, Verticies,snippet,[],[],'no');
        
        aaas=aaas+1;
    end
    
    
    % RevPet Primary Melts
    aaas=1;
    for zz = UniqueLabels_Rows'
        i =find(ismember(CombinedLabel, CombinedLabel(zz)));
        p=ismember(i,highlight_ReverseFC);
        i(p==0)=[];
        
        snippet =  ['''',runningMarkersPrimary{zz},''',''linewidth'',1,''MarkerSize'',',runningSizesPrimary{zz},',''MarkerEdgeColor'',', runningColors_LINE{zz},',''MarkerFaceColor'',', runningColors_FILL{zz}];
        GenericTernaryPlotter5(BestFitPrimaryTormey_n(i,:), IndiciesHere, Verticies,snippet,[],[],'no');
        
        aaas=aaas+1;
    end
    
    
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    % plots axis
    for llkm=highlight_ReverseFC(1)
        snippet =  ['''',runningMarkers{llkm},''',''linewidth'',1,''MarkerSize'',1,''MarkerEdgeColor'',''k'',''MarkerFaceColor'',', runningColors_FILL{llkm},''',''HandleVisibility'',''off'''];
        %GenericTernaryPlotter5(runningTormey((llkm),:), IndiciesHere, Verticies,snippet,[],[],'yes');
        GenericTernaryPlotter5([0 1 0 0], IndiciesHere, Verticies,snippet,[],[],'yesContour');
    end
    
    
    
    
    axis equal
    alabel = Verticies{IndiciesHere(1)};
    blabel = Verticies{IndiciesHere(2)};
    clabel = Verticies{IndiciesHere(3)};
    %Top Corner
    text(0.5,sqrt(3)/2, {alabel},'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',12,'FontName','Times New Roman');
    %Left Corner
    text(0, 0, {blabel},'VerticalAlignment','top','HorizontalAlignment','right','FontSize',12,'FontName','Times New Roman');
    %Right Corner
    text(1,0,{clabel},'HorizontalAlignment','left','VerticalAlignment','top','FontSize',12,'FontName','Times New Roman');
    
    axis(axishere(numTernaryPlots,:))
    



    if numTernaryPlots == 1
%         l=legend({'Gale et al. 2013 MORB glasses','Reverse FC Path',UniqueLabels{:},UniqueLabels_Primary{:}...
%             ['\it{Small: OPALM er >',num2str(round(GoodDataError*100,2,'significant')),'%}']},...
%             'FontName','Times New Roman','FontSize',12);
 % ['\it{Small: OPALM er >',num2str(round(GoodDataError*100,2,'significant')),'%}']         
            l=legend({'Gale et al. 2013 MORB glasses','Reverse FC Path',Unique_BAD_labels_leg{:}, UniqueLabels_leg{:},UniqueLabels_Primary{:}},...
            'FontName','Times New Roman','FontSize',12);    
        
    end
    
    set(gcf,'Position', [106 27 1075 902])
    set(gca,'FontName','Times New Roman')
    
end

set(l, 'Position',[0.2449 0.4982 0.1062 0.3753],'color','none')
set(gcf,'Position',[88 470 1525 473])
set(gcf,'name','Ternaries')

% subaxis(1,3,1)
% %axis([0.1264    0.8446    0.2433    0.8652]);
% axis([ 0.2902    0.4842    0.4961    0.6640]);
%
% IndiciesHere = Ternaries2plot(1,:);
% alabel = Verticies{IndiciesHere(1)};
% blabel = Verticies{IndiciesHere(2)};
% clabel = Verticies{IndiciesHere(3)};
%  a=axis;
%  %Top Corner
% text(a(4)./sqrt(3),a(4), {alabel},'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',12,'FontName','Times New Roman');
%
% %Left Corner
% text(a(1), a(3), {blabel},'VerticalAlignment','top','HorizontalAlignment','right','FontSize',12,'FontName','Times New Roman');
% %Right Corner
% text(a(2),a(3),{clabel},'HorizontalAlignment','left','VerticalAlignment','top','FontSize',12,'FontName','Times New Roman');
%
%
% subaxis(1,3,2)
% %axis([ 0.0672    0.7769    0.1380    0.7526]);
% axis([0.1399    0.5118    0.1897    0.5117]);
%
% IndiciesHere = Ternaries2plot(2,:);
% alabel = Verticies{IndiciesHere(1)};
% blabel = Verticies{IndiciesHere(2)};
% clabel = Verticies{IndiciesHere(3)};
% a=axis;
% %Left Corner
% text(a(1), a(3), {blabel},'VerticalAlignment','top','HorizontalAlignment','right','FontSize',12,'FontName','Times New Roman');
% %Right Corner
% text(a(2),a(3),{clabel},'HorizontalAlignment','left','VerticalAlignment','top','FontSize',12,'FontName','Times New Roman');
% %Top Corner
% text(a(4)./sqrt(3),a(4), {alabel},'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',12,'FontName','Times New Roman');
%
%
%
%
% subaxis(1,3,3)
% IndiciesHere = Ternaries2plot(3,:);
% alabel = Verticies{IndiciesHere(1)};
% blabel = Verticies{IndiciesHere(2)};
% clabel = Verticies{IndiciesHere(3)};
%
% axis([0.2281    0.4979    0.0821    0.3158]);
% % axis([-0.0141    0.5875   -0.0247    0.4963]);
% a=axis;
% %Left Corner
% text(a(1), a(3), {blabel},'VerticalAlignment','top','HorizontalAlignment','right','FontSize',12,'FontName','Times New Roman');
%
% %Right Corner
% text(a(2),a(3),{clabel},'HorizontalAlignment','left','VerticalAlignment','top','FontSize',12,'FontName','Times New Roman');
% %Top Corner
% text(a(1)+a(4)./sqrt(3),a(4), {alabel},'HorizontalAlignment','center','VerticalAlignment','bottom','FontSize',12,'FontName','Times New Roman');
%



%% Generic Harker Diagrams
figure(20)
close
figure(20)
set(gcf, 'Units', 'Inches', 'Position',  [0.1389 0.9167 17.2222 8.8611], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
hold on


% Choose from:
% {'T(C)','P(kbars)','SiO2','TiO2','Al2O3','Cr2O3','FeO','MnO','MgO','CaO','Na2O','K2O','P2O5','NiO','H2O','total',...
%     'Mg#','NaK#','CaO/Al2O3 wt%','CaO/Al2O3 moles','Na2O/FeO','K2O/TiO2','1-Mg#',...
%     'Qtz','Plag','Oliv','Cpx','Ox','Or','Ap'}


numcolumns = 4;
numrows = 2;

IndiciesFig1_X = {...
    'CaO'  'CaO'  'CaO'  'TiO2'; ...
    'FeO'  'MgO' 'TiO2' 'Mg#'; ...
    };


IndiciesFig1_Y = {...
    'SiO2'  'MgO' 'Al2O3' 'K2O'; ...
    'Na2O' 'K2O' 'Mg#' 'NaK#'; ...
    };


%zoomed out
axesMajors = [...
    0 15 42 61;
    0 15 0 21 ;
    0 15 8 22 ;
    .5 5 .01 20
    4 20 0 5;
    0 20 0.01 10;
    0 5 .25 .75 ;
    .08 .9 .08 .7];

%zoomed in
axesMajors = [...
    6 14.5 47 61;
    6 14.5 0 14 ;
    6 14.5 10 19 ;
    .5 4.1 .01 3
    4 20 1.5 5;
    0 12 0.01 3;
    0 5 .25 .75 ;
    .08 .9 .08 .5];




IndiciesX = upper(IndiciesFig1_X);
IndiciesY = upper(IndiciesFig1_Y);
subfigcount =1;
for kk = 1:size(IndiciesX,1)
    for n = 1:size(IndiciesX,2)
        
        
        yvalue =  IndiciesY(kk,n);
        xvalue = IndiciesX(kk,n);
        
        [a,xvalue]=(ismember(xvalue, upper(Elements)));
        [a,yvalue]=(ismember(yvalue, upper(Elements)));
        
        
        %for 3x3
        subaxis(numrows,numcolumns,subfigcount,'Spacing',.03,'Margin',0.015,'MarginLeft',0.017,'Padding',0.00,'PaddingLeft',0.02,'PaddingBottom',.033)
        % %for 1:2
        % set(gcf, 'Units', 'Inches', 'Position',   [20.3194 6.0278 10.9722 6.3194],  'PaperSize', [8.5, 11])
        % subaxis(numrows,numcolumns,n,'Spacing',.08,'Margin',0.02,'MarginBottom',.06,'MarginLeft',.1,'Padding',0,'PaddingLeft',0.02)
        
        
        %axis(axesMajors(subfigcount,:))
        hold on
        axis square
        
        
        h1=plot(Galeetal2013_Glasses(:,xvalue), Galeetal2013_Glasses(:,yvalue) ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',2,'MarkerFaceColor',grey8);% GALA14 517	Galapagos
        plot(Galeetal2013_Avg(:,xvalue), Galeetal2013_Avg(:,yvalue) ,'xk','Color',grey4,'LineWidth',3,'MarkerSize',15,'MarkerFaceColor',grey4);% GALA14 517	Galapagos
        
        
        % Paths of reverse FC
        a=1;
        for llkm=highlight_ReverseFC
             AllCorrectedLiquids_temp = BestFit_FCPath_n{llkm};
            if a==1
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','on');

            else
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','off');
            end

            a=a+1;
        end
        
        

    
    % Original Data that RevPet does not find a Valid Solution
    aaas=1;
    for zztemp = 1:size(Unique_BAD_labels_Rows,1)
        Unique_BAD_label_here = Unique_BAD_labels{zztemp};
        i =find(ismember(CombinedLabel, Unique_BAD_label_here));
        p=ismember(i,BAD_ReverseFC);
        i(p==0)=[];
        
        zz = i(1);
        plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
            runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
        
        
        aaas=aaas+1;
    end
    
    
        % Original Data with Valid Solutions 
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');

            aaas=aaas+1;
        end
        
        
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_n(i,xvalue),BestFitPrimary_n(i,yvalue),runningMarkersPrimary{zz},'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
       
        
        [Xtext,Ytext,FigureTitle] = oxideLabel(xvalue,yvalue);
        xlabel([Xtext,' (wt%)'])
        ylabel([Ytext,' (wt%)'])
        
        
        box on
        %grid on
        %set(gca, 'Xtick',1500:10:1540)
        set(gca,'ticklength',2*[0.0200    0.0500])
        set(gca,'fontsize', 14,'LineWidth',.7, 'FontName','Times New Roman')
        
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        numIncre = (ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1))/increment;
        ax.YAxis.MinorTickValues = ax.XAxis.TickValues(1)-numIncre.*increment:increment:ax.XAxis.TickValues(end)+numIncre.*increment;
        
        
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
        if strcmp(Ytext,'Na_2O')
            increment = .5;
        end
        numIncre = (ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1))/increment;
        ax.YAxis.MinorTickValues = ax.YAxis.TickValues(1)-numIncre.*increment:increment:ax.YAxis.TickValues(end)+numIncre.*increment;
        
        
        if strcmp(Ytext,'K_2O')
            set(gca,'yscale','log')
            ax.YAxis.MinorTickValues=  [.02:.01:.1 .2:.1:1 2:1:10];
        end
        
        
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        FigureTitle = 'MajorElement';
        set(gcf,'name',regexprep(FigureTitle,'\_*',''))
        uistack(h1,'bottom')
        
        
        subfigcount = subfigcount+1;
    end
end



legendgoodies = {'Gale et al. 2013 Glasses','Average E-, N-, and D-MORB'};
hlegend = legend([legendgoodies 'Reverse FC Path' Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');
%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])




%% MgO-versus oxides
figure(21)
close
figure(21)
set(gcf, 'Units', 'Inches', 'Position', [4.3472 0.1111 8.9861 12.8889], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
hold on


% Choose from:
% {'T(C)','P(kbars)','SiO2','TiO2','Al2O3','Cr2O3','FeO','MnO','MgO','CaO','Na2O','K2O','P2O5','NiO','H2O','total',...
%     'Mg#','NaK#','CaO/Al2O3 wt%','CaO/Al2O3 moles','Na2O/FeO','K2O/TiO2','1-Mg#',...
%     'Qtz','Plag','Oliv','Cpx','Ox','Or','Ap'}


numcolumns = 2;
numrows = 3;

IndiciesFig1_X = {...
    'MgO'  'MgO'; ...
    'MgO'  'MgO'; ...
    'MgO'  'MgO'; ...
    };


IndiciesFig1_Y = {...
    'SiO2'  'FeO'; ...
    'Al2O3' 'TiO2'; ...
    'K2O' 'CaO'; ...
    };



axesMajors = [0 13 45 75; 0 13 7 14 ;...
    0 13 13 19 ; 0 13 0.4 2.5
    0 13 0.06 2; 0 13 6 13];




IndiciesX = upper(IndiciesFig1_X);
IndiciesY = upper(IndiciesFig1_Y);
subfigcount =1;
for kk = 1:size(IndiciesX,1)
    for n = 1:size(IndiciesX,2)
        
        
        yvalue =  IndiciesY(kk,n);
        xvalue = IndiciesX(kk,n);
        
        [a,xvalue]=(ismember(xvalue, upper(Elements)));
        [a,yvalue]=(ismember(yvalue, upper(Elements)));
        
        
        %for 3x3
        subaxis(numrows,numcolumns,subfigcount,'Spacing',.02,'Margin',0.015,'Padding',0.00,'PaddingLeft',0.04,'PaddingBottom',.033)
        
        % %for 1:2
        % set(gcf, 'Units', 'Inches', 'Position',   [20.3194 6.0278 10.9722 6.3194],  'PaperSize', [8.5, 11])
        % subaxis(numrows,numcolumns,n,'Spacing',.08,'Margin',0.02,'MarginBottom',.06,'MarginLeft',.1,'Padding',0,'PaddingLeft',0.02)
        
        
        %axis(axesMajors(subfigcount,:))
        hold on
        axis square
        
        
        h1=plot(Galeetal2013_Glasses(:,xvalue), Galeetal2013_Glasses(:,yvalue) ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',2,'MarkerFaceColor',grey8);% GALA14 517	Galapagos
        plot(Galeetal2013_Avg(:,xvalue), Galeetal2013_Avg(:,yvalue) ,'xk','Color',grey4,'LineWidth',3,'MarkerSize',15,'MarkerFaceColor',grey4);% GALA14 517	Galapagos
        
        
%         for zz = UniqueLabels_Rows'
%             i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%             
%             plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
%                 runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
%             
%         end
        
        
        % Paths of reverse FC
        a=1;
        for llkm=highlight_ReverseFC
             AllCorrectedLiquids_temp = BestFit_FCPath_n{llkm};
            if a==1
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','on');

            else
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','off');
            end

            a=a+1;
        end
        
        
        % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            
            aaas=aaas+1;
        end
        
        
        
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            
         
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
           
            aaas=aaas+1;
        end
        
        
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_n(i,xvalue),BestFitPrimary_n(i,yvalue),runningMarkersPrimary{zz},'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
        
        
        
        
        [Xtext,Ytext,FigureTitle] = oxideLabel(xvalue,yvalue);
        xlabel([Xtext,' (wt%)'])
        ylabel([Ytext,' (wt%)'])
        
        
        box on
        %grid on
        %set(gca, 'Xtick',1500:10:1540)
        set(gca,'ticklength',2*[0.0200    0.0500])
        set(gca,'fontsize', 15,'LineWidth',.7,'FontName','Times New Roman')
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
        if strcmp(Ytext,'Na_2O')
            increment = .5;
        end
        
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
        %
        if strcmp(Ytext,'K_2O')
            set(gca,'yscale','log')
            ax.YAxis.MinorTickValues=  [.02:.01:.1 .2:.1:1 2:1:10];
        end
        
        
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        %set(gcf,'position',[2.7917 3.4306 9.7361 9.3750])
        % OPTION: uncomment for single (non-subplot plots)
        % legend(TabNames,'location','best')
        % title(FigureTitle)
        % FigureTitle = regexprep(FigureTitle,'\_*','');
        % FigureTitle = regexprep(FigureTitle,'\/','-');
        FigureTitle = 'Major_Element_MgO';
        set(gcf,'name',regexprep(FigureTitle,'\_*',''))
        %uistack(h2,'bottom')
        uistack(h1,'bottom')
        
        
        subfigcount = subfigcount+1;
    end
end


legendgoodies = {'Gale et al. 2013 Glasses','Average E-, N-, and D-MORB'};
hlegend = legend([legendgoodies 'Reverse FC Path' Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');

%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])




%% Money Figure (FeO-MgO and FeO-Na2O)
figure(22)
close
figure(22)
set(gcf, 'Units', 'Inches', 'Position',   [4.3472 6.0556 13.5278 6.9444], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
hold on


% Choose from:
% {'T(C)','P(kbars)','SiO2','TiO2','Al2O3','Cr2O3','FeO','MnO','MgO','CaO','Na2O','K2O','P2O5','NiO','H2O','total',...
%     'Mg#','NaK#','CaO/Al2O3 wt%','CaO/Al2O3 moles','Na2O/FeO','K2O/TiO2','1-Mg#',...
%     'Qtz','Plag','Oliv','Cpx','Ox','Or','Ap'}


numcolumns = 2;
numrows = 1;

IndiciesFig1_X = {...
    'FeO'  'FeO'; ...
    };


IndiciesFig1_Y = {...
    'MgO'  'Na2O'; ...
    };


axesMajors = [0 25 0 40; 0 25 0 10];





IndiciesX = upper(IndiciesFig1_X);
IndiciesY = upper(IndiciesFig1_Y);
subfigcount =1;
for kk = 1:size(IndiciesX,1)
    for n = 1:size(IndiciesX,2)
        
        
        yvalue =  IndiciesY(kk,n);
        xvalue = IndiciesX(kk,n);
        
        [a,xvalue]=(ismember(xvalue, upper(Elements)));
        [a,yvalue]=(ismember(yvalue, upper(Elements)));
        
        
        %for 3x3
        subaxis(numrows,numcolumns,subfigcount,'Spacing',.02,'Margin',0.015,'Padding',0.00,'PaddingLeft',0.04,'PaddingBottom',.033)
        
        % %for 1:2
        % set(gcf, 'Units', 'Inches', 'Position',   [20.3194 6.0278 10.9722 6.3194],  'PaperSize', [8.5, 11])
        % subaxis(numrows,numcolumns,n,'Spacing',.08,'Margin',0.02,'MarginBottom',.06,'MarginLeft',.1,'Padding',0,'PaddingLeft',0.02)
        
        
        %axis(axesMajors(subfigcount,:))
        hold on
        axis square
        
        
        h1=plot(Galeetal2013_Glasses(:,xvalue), Galeetal2013_Glasses(:,yvalue) ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',2,'MarkerFaceColor',grey8);% GALA14 517	Galapagos
        plot(Galeetal2013_Avg(:,xvalue), Galeetal2013_Avg(:,yvalue) ,'xk','Color',grey4,'LineWidth',3,'MarkerSize',15,'MarkerFaceColor',grey4);% GALA14 517	Galapagos
        
%         
%         for zz = UniqueLabels_Rows'
%             i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%             
%             plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
%                 runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
%             
%         end
        
        
               % Paths of reverse FC
        a=1;
        for llkm=highlight_ReverseFC
             AllCorrectedLiquids_temp = BestFit_FCPath_n{llkm};
            if a==1
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','on');

            else
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','off');
            end

            a=a+1;
        end
        
          % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            
            aaas=aaas+1;
        end
        
        
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
   
            aaas=aaas+1;
        end
        
        
         aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_n(i,xvalue),BestFitPrimary_n(i,yvalue),runningMarkersPrimary{zz},'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
        
        
        
        [Xtext,Ytext,FigureTitle] = oxideLabel(xvalue,yvalue);
        xlabel([Xtext,' (wt%)'])
        ylabel([Ytext,' (wt%)'])
        
        
        box on
        %grid on
        %
        set(gca,'ticklength',1.5*[0.0200    0.0500])
        set(gca,'fontsize', 17,'LineWidth',.7,'FontName','Times New Roman')
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        if strcmp(Xtext,'FeO')
            %increment = 0.5;
            set(gca, 'Xtick',0:1:100)
        end
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        
        if strcmp(Ytext,'MgO')
            %increment = 0.5;
            set(gca, 'Ytick',0:2:100)
        end
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
        
        
        
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
        %
        if strcmp(Ytext,'K_2O')
            set(gca,'yscale','log')
            ax.YAxis.MinorTickValues=  [.02:.01:.1 .2:.1:1 2:1:10];
        end
        
        
        
        
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        %set(gcf,'position',[2.7917 3.4306 9.7361 9.3750])
        % OPTION: uncomment for single (non-subplot plots)
        % legend(TabNames,'location','best')
        % title(FigureTitle)
        % FigureTitle = regexprep(FigureTitle,'\_*','');
        % FigureTitle = regexprep(FigureTitle,'\/','-');
        FigureTitle = 'Major_Element_2Panel';
        set(gcf,'name',regexprep(FigureTitle,'\_*',''))
        %uistack(h2,'bottom')
        uistack(h1,'bottom')
        
        
        subfigcount = subfigcount+1;
    end
end

legendgoodies = {'Gale et al. 2013 Glasses','Average E-, N-, and D-MORB'};
hlegend = legend([legendgoodies 'Reverse FC Path' Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');
%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])



%% 4 panel with money figures
figure(23)
close
figure(23)
hold on

specialMarkersSizes = [ones(1,9)*8 13 9 13];
DataLabelsTrace= {'SiO2'	'TiO2'	'Al2O3'	'Cr2O3'	'FeO'	'MgO'	'CaO'	'K2O'	'Na2O' 'total'  'Mg#' 'NaK#' 'CaO/Al2O3wt'	 'Rb'	'Sr'	'Cs'	'Ba'	'Pb'	'Ti'	'Y'	'Zr'	'Nb'	'Hf'	'Ta'	'Th'	'U'	'La'	'Ce'	'Pr'	'Nd'	'Sm'	'Eu'	'Gd'	'Tb'	'Dy'	'Ho'	'Er'	'Yb'	'Lu'};
targetStrings_Majors = {'Temp','Pressure','SiO2','TiO2','Al2O3', 'Cr2O3','FeO','MnO','MgO','CaO','Na2O','K2O','P2O5','NiO','H2O','total','Mg#','NaK#','CaO/Al2O3wt','CaO/Al2O3moles' 'Na2O/FeO'};

DataLabels = {'SiO2'	'TiO2'	'Al2O3'	'Cr2O3'	'FeO'	'MgO'	'CaO'	'K2O'	'Na2O'	'total'  'Mg#' 'NaK#' 'CaO/Al2O3wt' 'Na2O/FeO'};
% XXX important figure here!
%colorsjet = colormap(jet(size(test,2)));
clear xindice

ColorsHERE



IndiciesFig1_X = {...
    'FeO'  'FeO'; ...
    'Mg#'  'MgO'; ...
    };


IndiciesFig1_Y = {...
    'MgO'  'Na2O'; ...
    'Al2O3' 'CaO'; ...
    };


% axesMajors = [3 11 6 14; 3 11 48 55  ;...
%     3 11 .5 3 ; 3 11 15 21
%     3 11 7 12; 3 11 0.2 1.4];



axesMajors = [...
    5 20 0 40; 0 20 0 10;
    0.15 0.9 4 22; 0 35 0 17];


IndiciesX = upper(IndiciesFig1_X);
IndiciesY = upper(IndiciesFig1_Y);

numcolumns = 2;
numrows = 2;%round(size(Indicies,2)/numcolumns);
zzz=1;
for kk = 1:size(IndiciesX,1)
    for n = 1:size(IndiciesX,2)
        
        
        yvalueX =  IndiciesY(kk,n);
        xvalueY = IndiciesX(kk,n);
        
        [a,yvalue]=(ismember(yvalueX, upper(Elements)));
        [a,xvalue]=(ismember(xvalueY, upper(Elements)));
        
        
        set(gcf, 'Units', 'Inches', 'Position',  [7.8611 5.7500 7.4028 7.4861], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        subaxis(numrows,numcolumns,zzz,'Spacing',.03,'Margin',0.015,'Padding',0.00,'PaddingLeft',0.06,'PaddingBottom',.033)
        
        axis(axesMajors(zzz,:))
        hold on
        axis square
        
        
        
        MgNum_all = [0.5 .6 .7 .8];
        % MgNum_all = [0.2 .4:.1:.8];
        if kk==1 & n == 1
            %plotMgNumContours(xvalueY, yvalueX,axesMajors(zzz,[1 2]),axesMajors(zzz,[3 4]),MgNum_all)
            plotMgNumContours(xvalueY, yvalueX,axesMajors(zzz,[1 2]),[4 30],MgNum_all)
        end
        
        
        
        h1=plot(Galeetal2013_Glasses(:,xvalue), Galeetal2013_Glasses(:,yvalue) ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',2,'MarkerFaceColor',grey8);% GALA14 517	Galapagos
        plot(Galeetal2013_Avg(:,xvalue), Galeetal2013_Avg(:,yvalue),'xk','MarkerSize',15, 'LineWidth',3,'Color',grey4,'MarkerFaceColor',grey4);
        
        
%         aaas=1;
%         for zz = UniqueLabels_Rows'
%             i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%             
%             plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
%                 runningMarkers{zz},'LineWidth',0.5,'MarkerSize',str2num(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
%             aaas=aaas+1;
%         end
%         

        % Paths of reverse FC
        a=1;
        for llkm=highlight_ReverseFC
             AllCorrectedLiquids_temp = BestFit_FCPath_n{llkm};
            if a==1
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','on');

            else
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','off');
            end

            a=a+1;
        end
        
        
             % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            
            aaas=aaas+1;
        end
        
        
        % Erupted Basalts
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
        
            aaas=aaas+1;
        end
        
        % Best-fit Primary Basalts
          aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_n(i,xvalue),BestFitPrimary_n(i,yvalue),runningMarkersPrimary{zz},'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
        
        %Elements
        NeededYKG = {'SiO2', 'TiO2','Al2O3','Fe2O3','FeO','MgO','CaO','Na2O','K2O','P2O5','Cr2O3','MnO'} ;
        [ssdx,ElementIndicies4Target]=ismember(NeededYKG,Elements);
        noData = find(ElementIndicies4Target==0);
        ElementIndicies4Target(ElementIndicies4Target == 0) = max(ElementIndicies4Target);
        newlyOrderedElements = runningMajors(:,ElementIndicies4Target);
        newlyOrderedElements(:,noData)=[NaN];
        
        
        for sdf = [];%[1 10]
            [testfrac,testXXX,testYYY,plag,dd,CHANGEINDICIES] = YKG_wTRACE_Magnetite(newlyOrderedElements(sdf,:), 0.001, runningTrace(sdf,:),ones(size(runningTrace,2),7));
            [A,ElementIndicies4Target_YKG] = ismember(Elements, NeededYKG);
            testfracRightOrder = sumMgNumNorm_PT(testfrac,targetStrings_Majors,ElementIndicies4Target_YKG);
            
            f= [testfracRightOrder(:,9)];
            val = 5;
            tmp = abs(f-val);
            [idx idy] = min(tmp); %index of closest value
            closest = idy;
            plot(testfracRightOrder(1:closest,xvalue), testfracRightOrder(1:closest,yvalue),'-','Color','k','LineWidth',1,'handlevisibility','off')%'Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',markersize)
            plot(testfracRightOrder(CHANGEINDICIES,xvalue), testfracRightOrder(CHANGEINDICIES,yvalue),'o','Color','k','MarkerFaceColor','k','MarkerSize',eval(runningSizesBad{1}),'handlevisibility','off')%'Color',airsuperiorityblue,'MarkerFaceColor',airsuperiorityblue,'MarkerSize',markersize)
            
            
        end
        
        
        
        
        [Xtext,Ytext,FigureTitle] = oxideLabel(xvalue,yvalue);
        xlabel([Xtext,' (wt%)'])
        ylabel([Ytext,' (wt%)'])
        
        
        box on
        %grid on
        %set(gca, 'Xtick',1500:10:1540)
        set(gca,'ticklength',2*[0.0200    0.0500])
        set(gca,'fontsize', 13,'LineWidth',0.7, 'FontName','Times New Roman')
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
        if strcmp(Ytext,'Na_2O')
            increment = .5;
        end
        
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
        
        if strcmp(Ytext,'K_2O')
            set(gca,'yscale','log')
            ax.YAxis.MinorTickValues=  [.02:.01:.1 .2:.1:1 2:1:10];
        end
        
        
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        %set(gcf,'position',[2.7917 3.4306 9.7361 9.3750])
        % OPTION: uncomment for single (non-subplot plots)
        % legend(TabNames,'location','best')
        % title(FigureTitle)
        % FigureTitle = regexprep(FigureTitle,'\_*','');
        % FigureTitle = regexprep(FigureTitle,'\/','-');
        FigureTitle = 'Major_Element_4PanelFC';
        set(gcf,'name',regexprep(FigureTitle,'\_*',''))
        %uistack(h2,'bottom')
        uistack(h1,'bottom')
        
        if strcmp(xvalueY,'K2O')==1
            set(gca,'yscale','log')
        end
        zzz=zzz+1;
    end
end


legendgoodies = {'Gale et al. 2013 Glasses','Average E-, N-, and D-MORB'};
hlegend = legend([legendgoodies 'Reverse FC Path' Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');
%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])



%% Miyrashiro plot

figure(24)
close
figure(24)
set(gcf, 'Units', 'Inches', 'Position', [4.5833 4.0417 8.8333 8.5000], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
hold on



IndiciesFig1_X = {...
    'SiO2'
    };

IndiciesFig1_Y = {...
    'FeO'  'MgO'};


IndiciesX = [IndiciesFig1_X];
IndiciesY = [IndiciesFig1_Y];

%zoomed out
axesMajors = [47 66 0 7];
axesMajors = [45 70 0 23];
%zoomed in
% axesMajors = [4 15 44 54; 2 15 5 21 ;...
%   2 15 8 18 ; 0 4.5 .01 6
%     4 20 1 5; 0 21 0.01 6;
%     4 20 0 21 ; .2 .8 .1 .35];


IndiciesX = upper(IndiciesFig1_X);
IndiciesY = upper(IndiciesFig1_Y);

numcolumns = 1;
numrows = 1;
for kk = 1:size(IndiciesX,1)
    for n = 1:size(IndiciesX,2)
        yvalue1 =  IndiciesY(kk,n);
        yvalue2 =  IndiciesY(kk,n+1);
        xvalue = IndiciesX(kk,n);
        
        [a,xvalue]=(ismember(xvalue, upper(Elements)));
        [a,yvalue1]=(ismember(yvalue1, upper(Elements)));
        [a,yvalue2]=(ismember(yvalue2, upper(Elements)));
        
        
        
        
        
        axis(axesMajors(n,:))
        hold on
        axis square
        
        h1=plot(Galeetal2013_Glasses(:,xvalue), Galeetal2013_Glasses(:,yvalue1)./Galeetal2013_Glasses(:,yvalue2),'o','Color',grey8,'LineWidth',0.1,'MarkerSize',2,'MarkerFaceColor',grey8);% GALA14 517	Galapagos
        plot(Galeetal2013_Avg(:,xvalue), Galeetal2013_Avg(:,yvalue1)./Galeetal2013_Avg(:,yvalue2) ,'xk','Color',grey4,'LineWidth',3,'MarkerSize',15,'MarkerFaceColor',grey4);% GALA14 517	Galapagos
        
        
        
%         for zz = UniqueLabels_Rows'
%             i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%             
%             plot(runningMajors(i,xvalue),runningMajors(i,yvalue1)./runningMajors(i,yvalue2),...
%                 runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
%             
%         end
        


 % Paths of reverse FC
        a=1;
        for llkm=highlight_ReverseFC
             AllCorrectedLiquids_temp = BestFit_FCPath_n{llkm};
            if a==1
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue1)./AllCorrectedLiquids_temp(:,yvalue2),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','on');

            else
                plot(AllCorrectedLiquids_temp(:,xvalue),AllCorrectedLiquids_temp(:,yvalue1)./AllCorrectedLiquids_temp(:,yvalue2),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','off');
            end

            a=a+1;
        end
        
            % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue1)./runningMajors(i,yvalue2),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            
            aaas=aaas+1;
        end
        
        
        % Erupted Basalts
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));

            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningMajors(i,xvalue),runningMajors(i,yvalue1)./runningMajors(i,yvalue2),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
        
            aaas=aaas+1;
        end
        
        % Best-fit Primary Basalts
          aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_n(i,xvalue),BestFitPrimary_n(i,yvalue1)./BestFitPrimary_n(i,yvalue2),runningMarkersPrimary{zz},'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
        h3=refline(1/6.4, -42.8./6.4);
        set(h3,'color','k')
        
        [Xtext,Ytext,FigureTitle] = oxideLabel(xvalue,yvalue);
        xlabel([Xtext,' (wt%)'])
        ylabel([' FeO/MgO(wt%)'])
        
        
        box on
        %grid on
        %set(gca, 'Xtick',1500:10:1540)
        set(gca,'ticklength',2*[0.0200    0.0500])
        set(gca,'fontsize', 17,'LineWidth',.7,'FontName','Times New Roman')
        
        
        
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        FigureTitle = 'Miyashiro';
        set(gcf,'name',regexprep(FigureTitle,'\_*',''))
        uistack(h1,'bottom')
        if ismember(n,[4 6])
            set(gca,'yscale','log')
        end
        
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement(ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
        
        
        if strcmp(Ytext,'Na_2O')
            increment = .5;
        end
        
        %ax.YAxis.MinorTickValues = .5:increment:ax.YLim(2);
        
        if strcmp(Ytext,'K_2O')
            ax.YAxis.MinorTickValues=  [.02:.01:.1 .2:.1:1 2:1:10];
        end
        
        
    end
    
end



legendgoodies = {'Gale et al. 2013 Glasses','Average E-, N-, and D-MORB'};
hlegend = legend([legendgoodies 'Reverse FC Path' Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');

%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])

%%  Total-Alkali Silica Diagram
figure(25)
close
figure(25)
set(gcf, 'Units', 'Inches', 'Position', [4.5833 4.0417 8.8333 8.5000], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
hold on


IndiciesFig1_X = {...
    'SiO2'
    };

IndiciesFig1_Y = {...
    'Na2O'  'K2O'};

IndiciesX = [IndiciesFig1_X];
IndiciesY = [IndiciesFig1_Y];

%zoomed out
%axesMajors = [47 56 1.5 5];
axesMajors = [38 60 0 14];
%zoomed in
% axesMajors = [4 15 44 54; 2 15 5 21 ;...
%   2 15 8 18 ; 0 4.5 .01 6
%     4 20 1 5; 0 21 0.01 6;
%     4 20 0 21 ; .2 .8 .1 .35];


IndiciesX = upper(IndiciesFig1_X);
IndiciesY = upper(IndiciesFig1_Y);
numcolumns = 1;
numrows = 1;%round(size(Indicies,2)/numcolumns);
for kk = 1:size(IndiciesX,1)%size(StableIndex,2)
    for n = 1:size(IndiciesX,2)%size(Indicies,2)
        
        yvalue1 =  IndiciesY(kk,n);
        yvalue2 =  IndiciesY(kk,n+1);
        xvalue = IndiciesX(kk,n);
        
        [a,xvalue]=(ismember(xvalue, upper(Elements)));
        [a,yvalue1]=(ismember(yvalue1, upper(Elements)));
        [a,yvalue2]=(ismember(yvalue2, upper(Elements)));
        
        
        axis(axesMajors(n,:))
        hold on
        axis square
        
        
        
        h1=plot(Galeetal2013_Glasses(:,xvalue), nansum([Galeetal2013_Glasses(:,yvalue1) Galeetal2013_Glasses(:,yvalue2)]'),'o','Color',grey8,'LineWidth',0.1,'MarkerSize',2,'MarkerFaceColor',grey8);% GALA14 517	Galapagos
        plot(Galeetal2013_Avg(:,xvalue), nansum([Galeetal2013_Avg(:,yvalue1) Galeetal2013_Avg(:,yvalue2)]') ,'xk','Color',grey4,'LineWidth',3,'MarkerSize',15,'MarkerFaceColor',grey4);% GALA14 517	Galapagos
        
        
        
        
%         for zz = UniqueLabels_Rows'
%             i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%             plot(runningMajors(i,xvalue),nansum([runningMajors(i,yvalue1) runningMajors(i,yvalue2)]'),...
%                 runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
%         end
%         
        
        % Paths of reverse FC
        a=1;
        for llkm=highlight_ReverseFC
             AllCorrectedLiquids_temp = BestFit_FCPath_n{llkm};
            if a==1
                plot(AllCorrectedLiquids_temp(:,xvalue),nansum([AllCorrectedLiquids_temp(:,yvalue1) AllCorrectedLiquids_temp(:,yvalue2)]'),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','on');

            else
                plot(AllCorrectedLiquids_temp(:,xvalue),nansum([AllCorrectedLiquids_temp(:,yvalue1) AllCorrectedLiquids_temp(:,yvalue2)]'),'-','LineWidth',2,'Color',...
                    eval(runningColors_FILL{llkm}),'HandleVisibility','off');
            end

            a=a+1;
        end
        
        
        
              % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
               
            plot(runningMajors(i,xvalue),nansum([runningMajors(i,yvalue1) runningMajors(i,yvalue2)]'),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            
            aaas=aaas+1;
        end
        
        
        % Erupted Basalts
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningMajors(i,xvalue),nansum([runningMajors(i,yvalue1) runningMajors(i,yvalue2)]'),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
        
            aaas=aaas+1;
        end
        
        % Best-fit Primary Basalts
          aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_n(i,xvalue),nansum([BestFitPrimary_n(i,yvalue1) BestFitPrimary_n(i,yvalue2)]'),runningMarkersPrimary{zz},'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
        
        
        [Xtext,Ytext,FigureTitle] = oxideLabel(xvalue,yvalue);
        xlabel([Xtext,' (wt%)'])
        ylabel(['Na_2O+K_2O(wt%)'])
        
        
        box on
        %grid on
        set(gca, 'Xtick',38:2:80)
        set(gca,'ticklength',2*[0.0200    0.0500])
        set(gca,'fontsize', 17,'LineWidth',.7,'FontName','Times New Roman')
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        
        
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        
        
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
        if strcmp(Ytext,'Na_2O')
            increment = .5;
        end
        
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
        
        if strcmp(Ytext,'K_2O')
            ax.YAxis.MinorTickValues=  [.02:.01:.1 .2:.1:1 2:1:10];
        end
        
        
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        %set(gcf,'position',[2.7917 3.4306 9.7361 9.3750])
        % OPTION: uncomment for single (non-subplot plots)
        % legend(TabNames,'location','best')
        % title(FigureTitle)
        % FigureTitle = regexprep(FigureTitle,'\_*','');
        % FigureTitle = regexprep(FigureTitle,'\/','-');
        FigureTitle = 'TAS';
        set(gcf,'name',regexprep(FigureTitle,'\_*',''))
        %uistack(h2,'bottom')
        uistack(h1,'bottom')
        if ismember(n,[4 6])
            set(gca,'yscale','log')
        end
        
        
    end
    
    
end




legendgoodies = {'Gale et al. 2013 Glasses','Average E-, N-, and D-MORB'};
hlegend = legend([legendgoodies 'Reverse FC Path' Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');

%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])



%% Major Element - Trace Element
figure(30)
close
figure(30)
hold on

IndiciesFig1_X = {...
    'TiO2'  'MgO'; ...
    'MgO'  'MgO'; ...
    };

IndiciesFig1_Y = {...
    'Sr'  'Sr'; ...
    'Ba' 'Cr'; ...
    };


% axesMajors = [3 11 6 12; 3 11 48 55  ;...
%     3 11 .5 3 ; 3 11 15 21
%     3 11 7 12; 3 11 0.2 1.4];

IndiciesX = upper(IndiciesFig1_X);
IndiciesY = upper(IndiciesFig1_Y);

numcolumns = 2;
numrows = 2;%round(size(Indicies,2)/numcolumns);
zza=1;
for kk = 1:size(IndiciesX,1)%size(StableIndex,2)
    for n = 1:size(IndiciesX,2)%size(Indicies,2)
        
        
        ytext =  IndiciesY(kk,n);
        xtext= IndiciesX(kk,n);
        
        [a,xvalue]=(ismember(xtext, upper(Elements)));
        [a,yvalue]=(ismember(ytext, upper(targetStrings_Trace)));
        
        
        
        figure(30)
        %set(gcf, 'Units', 'Inches', 'Position',  [3.7917 2.8333 7.4028 9.7778], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        set(gcf, 'Units', 'Inches', 'Position',  [3.7917 0.2500 12.7083 12.3611], 'PaperUnits', 'Inches', 'PaperSize', [8.5, 11])
        
        
        subaxis(numrows,numcolumns,zza,'Spacing',.03,'Margin',0.015,'Padding',0.00,'PaddingLeft',0.04,'PaddingBottom',.04)
        
        %axis(axesMajors(zz,:))
        hold on
        axis square
        
        h1=plot(Galeetal2013_MajorTrace(:,xvalue), Galeetal2013_MajorTrace_Trace(:,yvalue) ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',2,'MarkerFaceColor',grey8);% GALA14 517	Galapagos
           
        
        plot(Galeetal2013_Avg(:,xvalue),Galeetal2013_Avg_Trace(:,yvalue),...
        'xk','Color',grey4,'LineWidth',3,'MarkerSize',15,'MarkerFaceColor',grey4);% GALA14 517	Galapagos
    
        
        
%         for zz = UniqueLabels_Rows'
%             i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%             
%             plot(runningMajors(i,xindice),runningTrace(i,iiiy),...
%                 runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
%             
%         end
        

        
      % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
              plot(runningMajors(i,xvalue),runningTrace(i,yvalue),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            aaas=aaas+1;
        end
        
        
        
        % Erupted Basalts
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            
          
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningMajors(i,xvalue),runningTrace(i,yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
        
            aaas=aaas+1;
        end
        
        
        % Best-fit Primary Basalts
          aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_n(i,xvalue),BestFitPrimary_Trace_n(i,yvalue),runningMarkersPrimary{zz},'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
        
        xlabel([xtext{:},' (wt%)']);
        ylabel([ytext{:},' (ppm)']);
        
        box on
        
        ax = gca;
        ax.XAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
        ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
        ax.YAxis.MinorTick = 'on';
        increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
        if strcmp(yvalue,'NA2O')
            increment = .5;
        end
        
        ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
        
        if strcmp(yvalue,'K2O')
            ax.YAxis.MinorTickValues=  [.02:.01:.1 .2:.1:1 2:1:10];
        end
        
        set(gca,'ticklength',2*[0.0200    0.0500])
        
        set(gca,'fontsize', 17,'LineWidth',.7,'FontName','Times New Roman')
        set(gca,'XColor', 'k')
        set(gca,'YColor', 'k')
        FigureTitle = 'Major_TraceElement';
        set(gcf,'name',regexprep(FigureTitle,'\_*',''))
        zza=zza+1;
        
    end
end


legendgoodies = {'Gale et al. 2013 Glasses'};
hlegend = legend([legendgoodies Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');

%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])


%% Trace Element Ratio v. Trace Element Ratio
figure(40)
close
figure(40)
hold on



masterratiosX ={'Sm/Yb' 'Th/U' 'Lu/Hf' 'Sm/Yb'};
masterratiosY ={'La/Sm' 'Zr/Yb' 'Th/U' 'Nb/Ta'};


axesMinors = [...
    4e-1 50 1e-1 13
    .2 4  .2 80
    .04 4   .2 3
    4e-1 50 .4 1.5 ];


masterratiosY ={'La/Sm' 'Eu/Eu*' 'Th/U' 'Nb/Ta'};
masterratiosX ={'Sm/Yb' 'La/Sm' 'Lu/Hf' 'Zr/Yb'};


numrows=2;
numcolumns=2;
set(gcf,'Position',[360 1 693 704])

nnn=1:4;


logon =1;
for n = nnn
    subaxis(numrows,numcolumns,n,'Spacing',.09,'Margin',0.02,'MarginBottom',.07,'MarginLeft',.08,'Padding',0,'PaddingLeft',0.02)
    
    %subaxis(numrows,numcolumns,n,'Spacing',.08,'Margin',0.02,'MarginBottom',.06,'MarginLeft',.06,'Padding',0,'PaddingLeft',0.02)
    %axis(axesMinors(n,:))
    hold on
    axis square
    
    
    Desired_X=masterratiosX{n};
    Desired_Y=masterratiosY{n};
    [a,xvalue] = ismember(Desired_X,RatioLabels);
    [b,yvalue] = ismember(Desired_Y,RatioLabels);
    
    
    
    
    h1=plot(Galeetal2013_MajorTrace_DataRatios(:,xvalue)./Sources_DataRatios(xvalue),Galeetal2013_MajorTrace_DataRatios(:,yvalue)./Sources_DataRatios(yvalue)...
        ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',4,'MarkerFaceColor',grey8);
    
    
    plot(Galeetal2013_Avg_DataRatios(:,xvalue)./Sources_DataRatios(xvalue),Galeetal2013_Avg_DataRatios(:,yvalue)./Sources_DataRatios(yvalue),...
        'xk','Color',grey4,'LineWidth',3,'MarkerSize',15,'MarkerFaceColor',grey4);% GALA14 517	Galapagos
    
    
%     for zz = UniqueLabels_Rows'
%         i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%         
%         plot(runningTraceRatios(i,xvalue)./Sources_DataRatios(xvalue),...
%             runningTraceRatios(i,yvalue)./Sources_DataRatios(yvalue),...
%             runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
%         
%     end
    
    

  % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
             plot(runningTraceRatios(i,xvalue)./Sources_DataRatios(xvalue),...
             runningTraceRatios(i,yvalue)./Sources_DataRatios(yvalue),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            aaas=aaas+1;
        end
        
        
        
    % Erupted Basalts
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningTraceRatios(i,xvalue)./Sources_DataRatios(xvalue),...
             runningTraceRatios(i,yvalue)./Sources_DataRatios(yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
        
            aaas=aaas+1;
        end
        
        
        % Best-fit Primary Basalts
          aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];

            plot(BestFitPrimary_TraceRatios_n(i,xvalue)./Sources_DataRatios(xvalue),...
                BestFitPrimary_TraceRatios_n(i,yvalue)/Sources_DataRatios(yvalue),...
                runningMarkersPrimary{zz}, 'LineWidth',1,...
                'MarkerSize',eval(runningSizesPrimary{zz}), 'MarkerFaceColor',eval(runningColors_FILL_Primary{zz}),'Color',eval(runningColors_LINE_Primary{zz}),'handlevisibility','on');
            
            aaas=aaas+1;
        end
        
    
        
    test1 = yline(1,'k-','handlevisibility','off'); %hline(1,'k-');
    test2=xline(1,'k-','handlevisibility','off');
    uistack(test1,'bottom')
    uistack(test2,'bottom')
    
    xlabel(sprintf('%s_{%s}',RatioLabels{xvalue}, Sources_Labels{2}));
    ylabel(sprintf('%s_{%s}',RatioLabels{yvalue}, Sources_Labels{2}));
    
    
    box on
    %grid on
    %set(gcf,'Position',[114 219 973 420])
    FigureTitle = sprintf('%s v. %s normalized to %s',RatioLabels{xvalue}, RatioLabels{yvalue}, Sources_Labels{2});
    FigureTitle = regexprep(FigureTitle,'\_*','');
    FigureTitle = regexprep(FigureTitle,'\/','-');
    
    FigureTitle = 'TraceRatios';
    set(gcf,'name',regexprep(FigureTitle,'\_*',''))
    %title(FigureTitle)
    set(gcf,'name',FigureTitle)
    axis square
    set(gca,'ticklength',2*[0.0200    0.0500])
    %set(gca,'ticklength',2*get(gca,'ticklength'))
    set(gca,'fontsize', 17,'LineWidth',.7,'FontName','Times New Roman')
    set(gca,'XColor', 'k')
    set(gca,'YColor', 'k')
    if logon ==1
        set(gca,'yscale','log')
        set(gca,'xscale','log')
    end
    
    %set(gca, 'Xtick',1500:10:1540)
    ax = gca;
    %  ax.XAxis.MinorTick = 'on';
    %     increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
    %  ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
    %set(ax,'XMinorTick','on')
    
    %set(gca, 'Ytick',[.4 .6 .8  1 1.2 1.4  1.6 1.8 2 2.2 3 ])
    
    
    if logon ==1
        set(gca,'yscale','log')
        set(gca,'xscale','log')
    end
    
    % if n == 2
    %
    %
    % ax.XAxis.MinorTickValues = [[axesMinors(n,1):.001:.01] [.02:.1:1] [2:1:axesMinors(n,2)]];
    % end
    
end


legendgoodies = {'Gale et al. 2013 Glasses' 'Avg E-, N-, D-MORB'};
% hlegend = legend([legendgoodies UniqueLabels'  UniqueLabels_Primary],'location','Best');
hlegend = legend([legendgoodies Unique_BAD_labels_leg{:} UniqueLabels_leg{:} UniqueLabels_Primary{:}],'location','Best');

%set(hlegend, 'Position',[0.8459 0.3280 0.1484 0.3229])

%% Isotope-Isotope plot

stableratio = 1;
%targetStrings_Isotopes = {'87Sr/86Sr'	'143Nd/144Nd'	'206Pb/204Pb'	'207Pb/204Pb'	'208Pb/204Pb'	'176Hf/177Hf' '230Th/238U' '226Ra/230Th'};

ZOO = [0.702817	0.5129	21.199	15.767	40.382	0.28285
    0.703870	0.512900	19.150	15.590	39.030	0.283
    0.703060	0.512900	19.946	15.765	39.247	0.283
    0.706745	0.512650	19.245	15.596	39.555	0.2827
    0.705010	0.512450	17.636	15.481	38.850	0.2826
    0.70227	0.5132	18.39	15.506	38.029	0.2832
    0.703697	0.51305	18.344	15.469	38.013	0.2831];
ZOO_LABELS = {'HIMUA'
    'HIMUB'
    'FOZOB2'
    'EM2B'
    'EM1B'
    'MORB'
    'TH'};



figure(50)
close


masterratioX = [3 3 1 2  3 1 ];
masterratioY = [1 4 6 6 5 2];
axesMinors = [17 22 .702 .708
    17 22 15.4 15.9
    0.702 .710 .2826 .2836
    .5122 .5134 .2826 .2836
    17 22 37 41
    .702 .710 .5122 .5134];


JACKSON_average = [0.703323	NaN	20.228	15.69	39.743	NaN
    0.702863	NaN	20.64	15.761	39.963	NaN
    0.703223	NaN	19.492	15.593	39.137	NaN
    0.703365	NaN	19.138	15.621	39.053	NaN
    0.703165	NaN	19.495	15.584	39.234	NaN
    0.70295	NaN	19.349	15.543	39.105	NaN
    0.703081	NaN	19.411	15.571	39.229	NaN
    0.703639	NaN	19.458	15.585	39.401	NaN
    0.703969	NaN	19.628	15.637	39.478	NaN
    0.705588	NaN	19.21	15.604	39.408	NaN
    0.704683	NaN	19.073	15.581	38.795	NaN
    0.703988	NaN	19.305	15.579	39.115	NaN
    0.703874	NaN	18.831	15.581	38.896	NaN
    0.705154	NaN	18.27	15.546	38.77	NaN
    0.70454	NaN	17.866	15.498	38.837	NaN
    0.703123	NaN	19.077	15.569	38.716	NaN
    0.703697	NaN	18.344	15.469	38.013	NaN
    0.703126	NaN	18.419	15.458	38.048	NaN
    0.702817	0.5129	21.199	15.767	40.382	NaN
    0.708443	NaN	19.149	15.623	39.491	NaN
    0.704642	NaN	17.826	15.496	38.855	NaN
    0.70227	NaN	18.39	15.506	38.029	NaN
    0.703870	0.512900	19.150	15.590	39.030	0.283
    0.706745	0.512650	19.245	15.596	39.555	NaN
    0.705010	0.512450	17.636	15.481	38.850	NaN
    0.702847	0.512900	21.622	15.819	40.546	0.283
    0.703060	0.512900	19.946	15.765	39.247	0.283];

Jackson_average_label={'HIMU'
    'HIMU'
    'HIMU'
    'HIMU'
    'HIMU'
    'HIMU'
    'HIMU'
    'HIMU'
    'EM2'
    'EM2'
    'EM2'
    'EM2'
    'EM2'
    'EM2'
    'EM1'
    'TH'
    'TH'
    'TH'
    'HIMU_A'
    'EM2_A'
    'EM1_A'
    'MORB'
    'HIMU_B'
    'EM2_B'
    'EM1_B'
    'HIMU_B2'
    'FOZO_B2'};


numrows=2;
numcolumns=3;
hold on

for n = 1:6
    figure(50)
    % set(gcf,'Position',[360 1 693 704])
    set(gcf,'Position',[87 0 1083 704])
    subaxis(numrows,numcolumns,n,'Spacing',.05,'Margin',0.02,'MarginBottom',.06,'MarginLeft',.04,'Padding',0,'PaddingLeft',0.02)
    %axis(axesMinors(n,:))
    hold on
    axis square
    
    
    xvalue = masterratioX(n);
    yvalue = masterratioY(n);
    
    %                 plot(DataRatios_GALEsheet(ultraslowspreadingMORBS,xvalue)./DataRatios_PUM(xvalue),...
    %                 DataRatios_GALEsheet(ultraslowspreadingMORBS,yvalue)./DataRatios_PUM(yvalue),...
    %                 'o', 'color','m', 'markerfacecolor','m','markersize',3)
    %
    %                 plot(DataRatios_GALEsheet(slowspreadingMORBS,xvalue)./DataRatios_PUM(xvalue),...
    %                 DataRatios_GALEsheet(slowspreadingMORBS,yvalue)./DataRatios_PUM(yvalue),...
    %                 'o', 'color','r', 'markerfacecolor','r','markersize',3)
    %
    %                       plot(DataRatios_GALEsheet(intermediatespreadingMORBS,xvalue)./DataRatios_PUM(xvalue),...
    %                 DataRatios_GALEsheet(intermediatespreadingMORBS,yvalue)./DataRatios_PUM(yvalue),...
    %                 'o', 'color','b', 'markerfacecolor','b','markersize',3)
    %
    %                 plot(DataRatios_GALEsheet(fastspreadingMORBS,xvalue)./DataRatios_PUM(xvalue),...
    %                 DataRatios_GALEsheet(fastspreadingMORBS,yvalue)./DataRatios_PUM(yvalue),...
    %                 'o', 'color','c', 'markerfacecolor','c','markersize',3)
    
    
    %                 plot(Galeetal2013_MajorTrace_Isotopes(emorbs1,xvalue),...
    %                 Galeetal2013_MajorTrace_Isotopes(emorbs1,yvalue),...
    %                 'o', 'color',grey5, 'markerfacecolor',grey5,'markersize',6)
    %
    %                 plot(Galeetal2013_MajorTrace_Isotopes(highSmYbMORBS,xvalue),...
    %                 Galeetal2013_MajorTrace_Isotopes(highSmYbMORBS,yvalue),...
    %                 'o', 'color','r', 'markerfacecolor','r','markersize',4)
    %
    
    
    
    %%% XXXX ZOO
    
    modelGalePCA5 = [0.70281	0.51313	18.47200	15.51300	38.25600	0.28335
        0.70256	0.51313	18.43100	15.49399	38.01600	0.28322
        0.70415	0.51274	18.96111	15.61131	39.05659	0.28302
        0.70291	0.51296	19.11800	15.60300	38.95100	0.28303
        0.70388	0.51274	19.39510	15.63039	39.35398	0.28300
        0.70401	0.51272	19.31123	15.63101	39.31623	0.28298
        0.70292	0.51298	19.10600	15.58900	38.91100	0.28312
        0.70268	0.51304	18.98312	15.55587	38.64329	0.28315
        0.70272	0.51301	19.06808	15.58535	38.82484	0.28308
        0.70269	0.51303	19.06376	15.55773	38.70530	0.28316
        0.70270	0.51303	19.03016	15.56766	38.72468	0.28312
        0.70258	0.51310	18.56838	15.52519	38.23540	0.28316];
    
    modelStandishPCA4=[0.70273	0.51305	18.56411	15.51207	38.19843	0.28310
        0.70259	0.51309	18.44453	15.49564	38.00739	0.28312
        0.70330	0.51289	19.06386	15.58070	38.99685	0.28302
        0.70299	0.51296	19.11800	15.60300	38.95100	0.28314
        0.70363	0.51279	19.35000	15.62000	39.45400	0.28297
        0.70354	0.51281	19.29900	15.61448	39.36348	0.28298
        0.70302	0.51298	19.10136	15.59288	38.91294	0.28334
        0.70256	0.51305	18.97161	15.59946	38.64668	0.28307
        0.70272	0.51301	19.07200	15.61000	38.81900	0.28308
        0.70249	0.51305	19.03800	15.61700	38.71300	0.28302
        0.70249	0.51305	19.03800	15.61700	38.71300	0.28302
        0.70273	0.51305	18.58479	15.51601	38.22474	0.28309];
    
    
    
    
    
    %    targetStrings_Isotopes = {'87Sr/86Sr'	'143Nd/144Nd'	'206Pb/204Pb'	'207Pb/204Pb'	'208Pb/204Pb'	'176Hf/177Hf'};
    %        text(JACKSON_average(:,xvalue),JACKSON_average(:,yvalue),Jackson_average_label)
    %
    %HIMUBand FOZO are sim, and still sim to HIMU
    ZOO = [0.702817	0.5129	21.199	15.767	40.382	0.28285
        0.706745	0.512650	19.245	15.596	39.555	0.2827
        0.703697	0.51305	18.344	15.469	38.013	0.2831
        0.703760 0.512865 19.479 15.650 39.086 0.282948];
    ZOO_LABELS = {'HIMU'
        'EM2'
        'TH'
        'Bouvet'};
    plot(ZOO(:,xvalue),ZOO(:,yvalue),'ko','MarkerFaceColor','k');
    text(ZOO(:,xvalue),ZOO(:,yvalue),ZOO_LABELS,'FontSize',10);
    
    
    
    %#4 works well, which is  1115. low Pb, low Sr, high Nd     'END0113-026-001'    'EPRR067'
    % testDMORB = (find(Galeetal2013_MajorTrace_Isotopes(:,3)<18&Galeetal2013_MajorTrace_Isotopes(:,1)<.7025));
    % Cluster1X=Galeetal2013_MajorTrace_Isotopes(testDMORB,:);
    % Cluster1X_labels= cellstr(num2str([1:size(testDMORB,1)]'));
    
    %#8 is the outlier, which is 2299 ('TRI0122-003-004'    'MARR081')
    % 12 might be more appropriate, which is 3451
    %GaleDataTextTraceOrdered(3451,:) MOA8801-022-013'    'SEIR066'
    %testDMORB = (find(Galeetal2013_MajorTrace_Isotopes(:,6)>.2834));
    
    %MORB must be the lowest Sr, which is not in the right place (pb is too
    %low). the closest pts dont have hf data
    %[a,testDMORB] = ((min(Galeetal2013_MajorTrace_Isotopes(:,1))));
    %[testDMORB] = find(Galeetal2013_MajorTrace_Isotopes(:,1)<.7024);
    %42 seems to be close to "MORB" #2954 END0063-002-005	MARR195
    %(1 is 760):NHOCHEP-031-004	EPRR026
    
    % newzooIndicies  = [1115 3451 2954 3158];
    % NEWZOO_LABELS={'LowPbSr','HighHf','MORB', 'EM1'};
    % NewZoo = Galeetal2013_MajorTrace_Isotopes(newzooIndicies,:);
    %          plot(NewZoo(:,xvalue),NewZoo(:,yvalue),'ko','MarkerFaceColor','b')
    %          text(NewZoo(:,xvalue),NewZoo(:,yvalue),NEWZOO_LABELS,'FontSize',10,'FontWeight','bold','Color','b')
    % 0.70229	0.51328	17.59400	15.40100	37.05700	0.28333
    % 0.70281	0.51326	17.76300	15.42800	37.59100	0.28344
    % 0.70237	0.51315	18.31600	15.51100	37.83600	0.28319
    
    newzooIndicies  = [1115 3451 2954 3158];
    NEWZOO_LABELS={'LowPbSr','HighHf','MORB', 'EM1'};
    NewZoo = Galeetal2013_MajorTrace_Isotopes(newzooIndicies,:);
    plot(NewZoo(:,xvalue),NewZoo(:,yvalue),'ko','MarkerFaceColor','k')
    text(NewZoo(:,xvalue),NewZoo(:,yvalue),NEWZOO_LABELS,'FontSize',10,'Color','k')
    
    %targetStrings_Isotopes = {'87Sr/86Sr'	'143Nd/144Nd'	'206Pb/204Pb'	'207Pb/204Pb'	'208Pb/204Pb'	'176Hf/177Hf'};
    %10 could work for PREMA but...
    [testDMORB] = find(Galeetal2013_MajorTrace_Isotopes(:,2)>.5128 & Galeetal2013_MajorTrace_Isotopes(:,2)<.51295 & Galeetal2013_MajorTrace_Isotopes(:,6)>.2831);
    
    
    % newzooIndicies  = [1115.00000	3133.00000	3158.00000	2853.00000	3236.00000	2299.00000];
    newzooIndicies  = [1115.00000	3133.00000	2853.00000	3236.00000	2299.00000];
    newzooIndicies_LABELS={'G1' 'G2' 'G4' 'G5' 'G6'};
    plot(Galeetal2013_MajorTrace_Isotopes(newzooIndicies,xvalue),Galeetal2013_MajorTrace_Isotopes(newzooIndicies,yvalue),'kp','MarkerFaceColor','b','MarkerSize',10)
    text(Galeetal2013_MajorTrace_Isotopes(newzooIndicies,xvalue),Galeetal2013_MajorTrace_Isotopes(newzooIndicies,yvalue),newzooIndicies_LABELS,'FontSize',10,'FontWeight','bold','Color','b')
    
    
    
%     for zz = UniqueLabels_Rows'
%         llkm =find(ismember(CombinedLabel, CombinedLabel(zz)));
%         
%         plot(runningIsotopes(llkm,xvalue),...
%             runningIsotopes(llkm,yvalue),...
%             runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}));
%         
%     end
    


  % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
            plot(runningIsotopes(i,xvalue),...
             runningIsotopes(i,yvalue),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            aaas=aaas+1;
        end
        
        
        
     % Erupted Basalts
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningIsotopes(i,xvalue),...
             runningIsotopes(i,yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
        
            aaas=aaas+1;
        end
        
      
        
    
    
    % modelStandishPCA4
    % modelGalePCA5
    % clear R_NRMSD_modelStandishPCA4
    % realData_TORMEY=CombinedStandish_Data_Isotopes(StandishUTh,1:6);
    %         for i=1:size(realData_TORMEY,1)
    %        % tempReal_NRMSD_Tormey = repmat(realData_TORMEY(i,:), length(realData_TORMEY), 1);
    %         R_NRMSD_modelStandishPCA4(i,:) = sqrt(nansum(((realData_TORMEY(i,:)-modelStandishPCA4(i,:))./realData_TORMEY(i,:)).^2, 2)./size(realData_TORMEY,2));
    %         end
    %
    %
    % clear R_NRMSD_Tormey_modelGalePCA5
    % realData_TORMEY=CombinedStandish_Data_Isotopes(StandishUTh,1:6);
    %         for i=1:size(realData_TORMEY,1)
    %        % tempReal_NRMSD_Tormey = repmat(realData_TORMEY(i,:), length(realData_TORMEY), 1);
    %         R_NRMSD_Tormey_modelGalePCA5(i,:) = sqrt(nansum(((realData_TORMEY(i,:)-modelGalePCA5(i,:))./realData_TORMEY(i,:)).^2, 2)./size(realData_TORMEY,2));
    %         end
    
    
    
    
    %#6 is a decent EM1, 3158: EWI9309-007-006	MARR244
    %[testDMORB] = find(Galeetal2013_MajorTrace_Isotopes(:,4)>15.55 & Galeetal2013_MajorTrace_Isotopes(:,3)<18);
    
    % Cluster1X=Galeetal2013_MajorTrace_Isotopes(testDMORB,:);
    % Cluster1X_labels= cellstr(num2str([1:size(testDMORB,1)]'));
    %
    %          plot(Cluster1X(:,xvalue),Cluster1X(:,yvalue),...
    %           'o','Color','r','LineWidth',0.1,'MarkerSize',6,'MarkerFaceColor','r');
    %           text(Cluster1X(:,xvalue),Cluster1X(:,yvalue),Cluster1X_labels,'FontSize',10,'FontWeight','bold','Color','k')
    
    
    %           [Cluster2,Cluster2X] = kmeans(Galeetal2013_MajorTrace_Isotopes(:,1:5),6);
    %          [Cluster1,Cluster1X] = kmeans(Galeetal2013_MajorTrace_Isotopes,6)
    %          Cluster1X_labels={'1' '2' '3' '4' '5' '6'};
    %          plot(Cluster1X(:,xvalue),Cluster1X(:,yvalue),'r^','MarkerSize',10,'MarkerFaceColor','r')
    %            text(Cluster1X(:,xvalue),Cluster1X(:,yvalue),Cluster1X_labels,'FontSize',10,'FontWeight','bold','Color','r')
    %
    %            Cluster2X(:,6) = NaN.*Cluster2X(:,5);
    %            Cluster2X_labels={'1' '2' '3' '4' '5' '6'};
    %            plot(Cluster2X(:,xvalue),Cluster2X(:,yvalue),'bs','MarkerSize',10,'MarkerFaceColor','b')
    %              text(Cluster2X(:,xvalue),Cluster2X(:,yvalue),Cluster2X_labels,'FontSize',10,'FontWeight','bold','Color','b')
    %
    
    
    
    
    
    h1=plot(Galeetal2013_MajorTrace_Isotopes(:,xvalue),Galeetal2013_MajorTrace_Isotopes(:,yvalue)...
        ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',6,'MarkerFaceColor',grey8);
    
    
    ylabel( regexprep(targetStrings_Isotopes{yvalue},'\_*','/'));
    xlabel( regexprep(targetStrings_Isotopes{xvalue},'\_*','/'));
    
    
    ax = gca;
    ax.XAxis.MinorTick = 'on';
    increment = findBestIncrement( ax.XLim(2)- ax.XLim(1),ax.XAxis.TickValues(2) -  ax.XAxis.TickValues(1));
    if (ax.XAxis.TickValues(2)- ax.XAxis.TickValues(1)) < increment
        increment = (ax.XAxis.TickValues(2)- ax.XAxis.TickValues(1))./2;
    end
    ax.XAxis.MinorTickValues = ax.XLim(1):increment:ax.XLim(2);
    
    
    ax.YAxis.MinorTick = 'on';
    increment = findBestIncrement( ax.YLim(2)- ax.YLim(1),ax.YAxis.TickValues(2) -  ax.YAxis.TickValues(1));
    if (ax.YAxis.TickValues(2)- ax.YAxis.TickValues(1)) < increment
        increment = (ax.YAxis.TickValues(2)- ax.YAxis.TickValues(1))./2;
    end
    ax.YAxis.MinorTickValues = ax.YLim(1):increment:ax.YLim(2);
    
    
    
    box on
    %grid on
    %set(gcf,'Position',[114 219 973 420])
    FigureTitle = 'Isotope_Figure';
    set(gcf,'name',regexprep(FigureTitle,'\_*',''))
    %title(FigureTitle)
    set(gcf,'name',FigureTitle)
    axis square
    set(gca,'ticklength',2*[0.0200    0.0500])
    %set(gca,'ticklength',2*get(gca,'ticklength'))
    set(gca,'fontsize', 13,'LineWidth',.7,'FontName','Times New Roman')
    set(gca,'XColor', 'k')
    set(gca,'YColor', 'k')
    
    set(gca,'yscale','log')
    set(gca,'xscale','log')
    
    
    uistack(h1,'bottom')
    
    set(gca,'yscale','log')
    set(gca,'xscale','log')
end




legendgoodies ={'GALE et al. 2013 MORB glasses' 'Zoo','Zoo','Zoo'};
hlegned = legend([legendgoodies  Unique_BAD_labels_leg{:} UniqueLabels_leg{:}],'location','Best');
set(hlegned,'Position',[0.8061 0.7969 0.1768 0.1641])



%% Trace Ratio v. Isotope


figure(60)
close


% options are:
% 87Sr/86Sr	143Nd/144Nd	206Pb/204Pb	207Pb/204Pb	208Pb/204Pb	176Hf/177Hf
masterratiosY ={'SR87_SR86' 'SR87_SR86' 'HF176_HF177' 'HF176_HF177'};
masterratiosX ={'Sm/Yb' 'La/Sm' 'Lu/Hf' 'Zr/Yb'};



numrows=2;
numcolumns=2;
hold on
for n = 1:4
    figure(60)
    %  set(gcf,'Position',[360 1 693 704])
    set(gcf,'Position', [360 1 741 704])
    
    subaxis(numrows,numcolumns,n,'Spacing',.08,'Margin',0.02,'MarginBottom',.06,'MarginLeft',.08,'Padding',0,'PaddingLeft',0.02)
    %axis(axesMinors(n,:))
    hold on
    axis square
    
    
    Desired_X=masterratiosX{n};
    Desired_Y=masterratiosY{n};
    [a,xvalue] = ismember(Desired_X,RatioLabels);
    [b,yvalue] = ismember(Desired_Y,targetStrings_Isotopes);
    
    
    h1=plot(Galeetal2013_MajorTrace_DataRatios(:,xvalue)./Sources_DataRatios(xvalue),Galeetal2013_MajorTrace_Isotopes(:,yvalue)...
        ,'o','Color',grey8,'LineWidth',0.1,'MarkerSize',6,'MarkerFaceColor',grey8);
    
%     for zz = UniqueLabels_Rows'
%         i =find(ismember(CombinedLabel, CombinedLabel(zz)));
%         
%         plot(runningTraceRatios(i,xvalue)./Sources_DataRatios(xvalue),...
%             runningIsotopes(i,yvalue),...
%             runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizes{zz}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}));
%     end
    
    
  % Original Data that RevPet does not find a Valid Solution
        aaas=1;
        for zztemp = 1:size(Unique_BAD_labels_Rows,1)
            Unique_BAD_label_here = Unique_BAD_labels{zztemp};
            i =find(ismember(CombinedLabel, Unique_BAD_label_here));
            p=ismember(i,BAD_ReverseFC);
            i(p==0)=[];
            
            zz = i(1);
            plot(runningTraceRatios(i,xvalue)./Sources_DataRatios(xvalue),...
             runningIsotopes(i,yvalue),...
                runningMarkers{zz},'LineWidth',0.5,'MarkerSize',eval(runningSizesBad{1}), 'Color',eval(runningColors_LINE{zz}),'MarkerFaceColor',eval(runningColors_FILL{zz}),'HandleVisibility','on')
            
            aaas=aaas+1;
        end
        
        
        
    
       % Erupted Basalts
        aaas=1;
        for zz = UniqueLabels_Rows'
            i =find(ismember(CombinedLabel, CombinedLabel(zz)));
            
          
            p=ismember(i,highlight_ReverseFC);
            i(p==0)=[];
            
            plot(runningTraceRatios(i,xvalue)./Sources_DataRatios(xvalue),...
             runningIsotopes(i,yvalue),...
                runningMarkers{zz},'LineWidth',1,'MarkerSize',eval(runningSizes{zz}), 'MarkerFaceColor',eval(runningColors_FILL{zz}),'Color',eval(runningColors_LINE{zz}),...
                'handlevisibility','on');
        
            aaas=aaas+1;
        end
        
        
        
    test2=xline(1,'k-','handlevisibility','off');
    
    uistack(test2,'bottom')
    
    
    
    xlabel(sprintf('%s normalized to %s',Desired_X, Sources_Labels{2}));
    ylabel( regexprep(Desired_Y,'\_*','/'));
    
    
    
    box on
    %grid on
    %set(gcf,'Position',[114 219 973 420])
    ax=gca;
    ax.YAxis.MinorTick = 'on';
    ax.XAxis.MinorTick = 'on';
    FigureTitle = 'TraceRatio_Isotope';
    set(gcf,'name',regexprep(FigureTitle,'\_*',''))
    %title(FigureTitle)
    set(gcf,'name',FigureTitle)
    axis square
    set(gca,'ticklength',2*[0.0200    0.0500])
    %set(gca,'ticklength',2*get(gca,'ticklength'))
    set(gca,'fontsize', 15,'LineWidth',.7,'FontName','Times New Roman')
    set(gca,'XColor', 'k')
    set(gca,'YColor', 'k')
    %  set(gca,'yscale','log')
    % set(gca,'xscale','log')
    
    
    uistack(h1,'bottom')
end




legendgoodies ={'GALE et al. 2013 MORB glasses'};
hlegned = legend([legendgoodies  Unique_BAD_labels_leg{:} UniqueLabels_leg{:}],'location','Best');
set(hlegned,'Position',[0.8061 0.7969 0.1768 0.1641])

return
%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%cab is a "close all but" function, so tell it the figuers you want to keep (once you have the axes the way you want) 
%and then run this section to save the figures e.g., 
cab 1 25
savefigsPDF(worksheetName2Save,subfolder)